In this brief example, I show how to use the GPS coordinates from the Demographic and Health Survey data and merge them to the ADM2 subnational geographic level for the country of Ethopia. Then I produce estimates using the DHS data for ADM 2 regions of the country.

This is possible by useing the GIS capacity of the sf package to spatially intersect the DHS points and the ADM 2 polygons.

library(sf)
## Warning: package 'sf' was built under R version 4.1.3
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(mapview)

Read in dhs points

ethpoly <- st_read(dsn = "C:/Users/ozd504//OneDrive - University of Texas at San Antonio//students/fikre/spatial_epi/ETH_adm2.shp")
## Reading layer `ETH_adm2' from data source 
##   `C:\Users\ozd504\OneDrive - University of Texas at San Antonio\students\fikre\spatial_epi\ETH_adm2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 72 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 32.99994 ymin: 3.322099 xmax: 47.98618 ymax: 14.89996
## Geodetic CRS:  WGS 84
ethpoly$struct <- 1:dim(ethpoly)[1]

plot(ethpoly["struct"])

Read in dhs sample locations and ADM 2 regions.

The adm2 shapefile can be found in the Diva GIS international data repository, or from the IPUMS International site below I use the ADM2 level of administrative geography.

These locations are not identified in the DHS, but by performing a spatial intersection, we can merge the DHS survey locations to the ADM 2 units

eth_dots<-st_read("C:/Users/ozd504//OneDrive - University of Texas at San Antonio//students//fikre/ethiopia_dhs/ETGE52FL/ETGE52FL.shp")
## Reading layer `ETGE52FL' from data source 
##   `C:\Users\ozd504\OneDrive - University of Texas at San Antonio\students\fikre\ethiopia_dhs\ETGE52FL\ETGE52FL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 535 features and 20 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 0 ymin: 0 xmax: 43.36409 ymax: 14.50379
## Geodetic CRS:  WGS 84
eth_dots <- eth_dots[eth_dots$LATNUM>0,]
eth_adm2<-st_read("C:/Users/ozd504//OneDrive - University of Texas at San Antonio//students//fikre/spatial_epi/ETH_adm2.shp")
## Reading layer `ETH_adm2' from data source 
##   `C:\Users\ozd504\OneDrive - University of Texas at San Antonio\students\fikre\spatial_epi\ETH_adm2.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 72 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 32.99994 ymin: 3.322099 xmax: 47.98618 ymax: 14.89996
## Geodetic CRS:  WGS 84
#merge dots to administrative data
eth_dots2000<-st_intersection(eth_dots, eth_adm2)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
mapview(eth_dots["DHSCLUST"])+mapview(eth_adm2["NAME_2"])

Read in DHS survey and recode stunting outcome

Next, I use the DHS survey data to estimate the prevalence of stunting in the ADM 2 regions.

library(haven)
dhs2000<-read_dta("C:/Users/ozd504//OneDrive - University of Texas at San Antonio//students//fikre/ethiopia_dhs/ETKR41DT/ETKR41FL.DTA")
dhs2000<-zap_labels(dhs2000)

library(car)
## Loading required package: carData
dhs2000$stunting<-ifelse(dhs2000$hw5/100<=-2&dhs2000$hw5/100!=-2,1,0)
#dhs2000$sex<-dhs2000$hc27

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
dhs2000<-dhs2000%>%
  mutate(wt = v005/1000000)%>%
  filter(complete.cases(stunting))%>%
  select(v001,stunting, wt, v021, v022)

Merge survey data to sample locations

dhs2000m<-merge(dhs2000, eth_dots2000, by.x="v001", by.y="DHSCLUST")

Create survey estimates for new regions after spatial intersection

library(survey)
## Loading required package: grid
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 4.1.3
## Loading required package: survival
## Warning: package 'survival' was built under R version 4.1.3
## 
## Attaching package: 'survey'
## The following object is masked from 'package:graphics':
## 
##     dotchart
options(survey.lonely.psu = "adjust")
des<-svydesign(ids = ~v021, strata = ~v022, weights = ~wt, data=dhs2000m)
names(dhs2000m)
##  [1] "v001"       "stunting"   "wt"         "v021"       "v022"      
##  [6] "DHSID"      "DHSCC"      "DHSYEAR"    "CCFIPS"     "ADM1FIPS"  
## [11] "ADM1FIPSNA" "ADM1SALBNA" "ADM1SALBCO" "ADM1DHS"    "ADM1NAME"  
## [16] "DHSREGCO"   "DHSREGNA"   "SOURCE"     "URBAN_RURA" "LATNUM"    
## [21] "LONGNUM"    "ALT_GPS"    "ALT_DEM"    "DATUM"      "ID_0"      
## [26] "ISO"        "NAME_0"     "ID_1"       "NAME_1"     "ID_2"      
## [31] "NAME_2"     "TYPE_2"     "ENGTYPE_2"  "NL_NAME_2"  "VARNAME_2" 
## [36] "geometry"
est.stunt <- svyby(~stunting, ~ID_2, des, FUN=svymean, na.rm=T)

head(est.stunt)

merge estimates to map and map stunting prevalence

library(tigris)
## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.
library(mapview)
mdat<- geo_join(ethpoly, est.stunt, "ID_2","ID_2")
## Warning: We recommend using the dplyr::*_join() family of functions instead.
## Warning: `group_by_()` was deprecated in dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
mapview(mdat["stunting"])
LS0tDQp0aXRsZTogIk1hcHBpbmcgREhTIFN1cnZleSBFc3RpbWF0ZXMgdG8gQURNIExldmVsIDIgR2VvZ3JhcGhpZXMiDQphdXRob3I6ICJDb3JleSBTcGFya3MgUGhEIg0KZGF0ZTogImByIGZvcm1hdChTeXMudGltZSgpLCAnJWQgJUIsICVZJylgIg0Kb3V0cHV0Og0KICAgaHRtbF9kb2N1bWVudDoNCiAgICBkZl9wcmludDogcGFnZWQNCiAgICBmaWdfaGVpZ2h0OiA3DQogICAgZmlnX3dpZHRoOiA3DQogICAgdG9jOiB5ZXMNCiAgICB0b2NfZmxvYXQ6IHllcw0KICAgIGNvZGVfZG93bmxvYWQ6IHRydWUNCi0tLQ0KDQpJbiB0aGlzIGJyaWVmIGV4YW1wbGUsIEkgc2hvdyBob3cgdG8gdXNlIHRoZSBHUFMgY29vcmRpbmF0ZXMgZnJvbSB0aGUgW0RlbW9ncmFwaGljIGFuZCBIZWFsdGggU3VydmV5IGRhdGFdKCkgYW5kIG1lcmdlIHRoZW0gdG8gdGhlIEFETTIgc3VibmF0aW9uYWwgZ2VvZ3JhcGhpYyBsZXZlbCBmb3IgdGhlIGNvdW50cnkgb2YgRXRob3BpYS4gVGhlbiBJIHByb2R1Y2UgZXN0aW1hdGVzIHVzaW5nIHRoZSBESFMgZGF0YSBmb3IgQURNIDIgcmVnaW9ucyBvZiB0aGUgY291bnRyeS4NCg0KVGhpcyBpcyBwb3NzaWJsZSBieSB1c2VpbmcgdGhlIEdJUyBjYXBhY2l0eSBvZiB0aGUgYHNmYCBwYWNrYWdlIHRvIHNwYXRpYWxseSBpbnRlcnNlY3QgdGhlIERIUyBwb2ludHMgYW5kIHRoZSBBRE0gMiBwb2x5Z29ucy4gDQoNCmBgYHtyfQ0KbGlicmFyeShzZikNCmxpYnJhcnkobWFwdmlldykNCmBgYA0KDQojIyBSZWFkIGluIGRocyBwb2ludHMNCmBgYHtyfQ0KZXRocG9seSA8LSBzdF9yZWFkKGRzbiA9ICJDOi9Vc2Vycy9vemQ1MDQvL09uZURyaXZlIC0gVW5pdmVyc2l0eSBvZiBUZXhhcyBhdCBTYW4gQW50b25pby8vc3R1ZGVudHMvZmlrcmUvc3BhdGlhbF9lcGkvRVRIX2FkbTIuc2hwIikNCg0KZXRocG9seSRzdHJ1Y3QgPC0gMTpkaW0oZXRocG9seSlbMV0NCg0KcGxvdChldGhwb2x5WyJzdHJ1Y3QiXSkNCg0KYGBgDQoNCg0KIyMgUmVhZCBpbiBkaHMgc2FtcGxlIGxvY2F0aW9ucyBhbmQgQURNIDIgcmVnaW9ucy4NCg0KVGhlIGFkbTIgc2hhcGVmaWxlIGNhbiBiZSBmb3VuZCBpbiB0aGUgW0RpdmEgR0lTIGludGVybmF0aW9uYWwgZGF0YSByZXBvc2l0b3J5XShodHRwczovL3d3dy5kaXZhLWdpcy5vcmcvZ2RhdGEpLCBvciBmcm9tIHRoZSBbSVBVTVMgSW50ZXJuYXRpb25hbCBzaXRlXShodHRwczovL2ludGVybmF0aW9uYWwuaXB1bXMub3JnL2ludGVybmF0aW9uYWwvZ2lzX2hhcm1vbml6ZWRfMm5kLnNodG1sKSBiZWxvdyBJIHVzZSB0aGUgQURNMiBsZXZlbCBvZiBhZG1pbmlzdHJhdGl2ZSBnZW9ncmFwaHkuIA0KDQpUaGVzZSBsb2NhdGlvbnMgYXJlIG5vdCBpZGVudGlmaWVkIGluIHRoZSBESFMsIGJ1dCBieSBwZXJmb3JtaW5nIGEgc3BhdGlhbCBpbnRlcnNlY3Rpb24sIHdlIGNhbiBtZXJnZSB0aGUgREhTIHN1cnZleSBsb2NhdGlvbnMgdG8gdGhlIEFETSAyIHVuaXRzDQoNCg0KDQpgYGB7cn0NCmV0aF9kb3RzPC1zdF9yZWFkKCJDOi9Vc2Vycy9vemQ1MDQvL09uZURyaXZlIC0gVW5pdmVyc2l0eSBvZiBUZXhhcyBhdCBTYW4gQW50b25pby8vc3R1ZGVudHMvL2Zpa3JlL2V0aGlvcGlhX2Rocy9FVEdFNTJGTC9FVEdFNTJGTC5zaHAiKQ0KZXRoX2RvdHMgPC0gZXRoX2RvdHNbZXRoX2RvdHMkTEFUTlVNPjAsXQ0KYGBgDQoNCmBgYHtyfQ0KZXRoX2FkbTI8LXN0X3JlYWQoIkM6L1VzZXJzL296ZDUwNC8vT25lRHJpdmUgLSBVbml2ZXJzaXR5IG9mIFRleGFzIGF0IFNhbiBBbnRvbmlvLy9zdHVkZW50cy8vZmlrcmUvc3BhdGlhbF9lcGkvRVRIX2FkbTIuc2hwIikNCmBgYA0KDQpgYGB7cn0NCiNtZXJnZSBkb3RzIHRvIGFkbWluaXN0cmF0aXZlIGRhdGENCmV0aF9kb3RzMjAwMDwtc3RfaW50ZXJzZWN0aW9uKGV0aF9kb3RzLCBldGhfYWRtMikNCg0KbWFwdmlldyhldGhfZG90c1siREhTQ0xVU1QiXSkrbWFwdmlldyhldGhfYWRtMlsiTkFNRV8yIl0pDQpgYGANCg0KIyMgUmVhZCBpbiBESFMgc3VydmV5IGFuZCByZWNvZGUgc3R1bnRpbmcgb3V0Y29tZQ0KDQpOZXh0LCBJIHVzZSB0aGUgREhTIHN1cnZleSBkYXRhIHRvIGVzdGltYXRlIHRoZSBwcmV2YWxlbmNlIG9mIHN0dW50aW5nIGluIHRoZSBBRE0gMiByZWdpb25zLg0KDQoNCmBgYHtyfQ0KbGlicmFyeShoYXZlbikNCmRoczIwMDA8LXJlYWRfZHRhKCJDOi9Vc2Vycy9vemQ1MDQvL09uZURyaXZlIC0gVW5pdmVyc2l0eSBvZiBUZXhhcyBhdCBTYW4gQW50b25pby8vc3R1ZGVudHMvL2Zpa3JlL2V0aGlvcGlhX2Rocy9FVEtSNDFEVC9FVEtSNDFGTC5EVEEiKQ0KZGhzMjAwMDwtemFwX2xhYmVscyhkaHMyMDAwKQ0KDQpsaWJyYXJ5KGNhcikNCmRoczIwMDAkc3R1bnRpbmc8LWlmZWxzZShkaHMyMDAwJGh3NS8xMDA8PS0yJmRoczIwMDAkaHc1LzEwMCE9LTIsMSwwKQ0KI2RoczIwMDAkc2V4PC1kaHMyMDAwJGhjMjcNCg0KbGlicmFyeShkcGx5cikNCmRoczIwMDA8LWRoczIwMDAlPiUNCiAgbXV0YXRlKHd0ID0gdjAwNS8xMDAwMDAwKSU+JQ0KICBmaWx0ZXIoY29tcGxldGUuY2FzZXMoc3R1bnRpbmcpKSU+JQ0KICBzZWxlY3QodjAwMSxzdHVudGluZywgd3QsIHYwMjEsIHYwMjIpDQpgYGANCg0KIyMgTWVyZ2Ugc3VydmV5IGRhdGEgdG8gc2FtcGxlIGxvY2F0aW9ucw0KDQpgYGB7cn0NCmRoczIwMDBtPC1tZXJnZShkaHMyMDAwLCBldGhfZG90czIwMDAsIGJ5Lng9InYwMDEiLCBieS55PSJESFNDTFVTVCIpDQoNCmBgYA0KDQoNCiMjIENyZWF0ZSBzdXJ2ZXkgZXN0aW1hdGVzIGZvciBuZXcgcmVnaW9ucyBhZnRlciBzcGF0aWFsIGludGVyc2VjdGlvbg0KDQpgYGB7cn0NCmxpYnJhcnkoc3VydmV5KQ0Kb3B0aW9ucyhzdXJ2ZXkubG9uZWx5LnBzdSA9ICJhZGp1c3QiKQ0KZGVzPC1zdnlkZXNpZ24oaWRzID0gfnYwMjEsIHN0cmF0YSA9IH52MDIyLCB3ZWlnaHRzID0gfnd0LCBkYXRhPWRoczIwMDBtKQ0KbmFtZXMoZGhzMjAwMG0pDQplc3Quc3R1bnQgPC0gc3Z5YnkofnN0dW50aW5nLCB+SURfMiwgZGVzLCBGVU49c3Z5bWVhbiwgbmEucm09VCkNCg0KaGVhZChlc3Quc3R1bnQpDQoNCmBgYA0KDQojIyBtZXJnZSBlc3RpbWF0ZXMgdG8gbWFwIGFuZCBtYXAgc3R1bnRpbmcgcHJldmFsZW5jZQ0KDQpgYGB7cn0NCmxpYnJhcnkodGlncmlzKQ0KbGlicmFyeShtYXB2aWV3KQ0KbWRhdDwtIGdlb19qb2luKGV0aHBvbHksIGVzdC5zdHVudCwgIklEXzIiLCJJRF8yIikNCm1hcHZpZXcobWRhdFsic3R1bnRpbmciXSkNCg0KDQpgYGA=